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Abstract 

A one - parameter dynamical system is associated to the mathe- 
matical problem governing the membrane excitability of a ventricular 
cardiomyocyte, according to the Luo-Rudy I model. An algorithm used 
to construct the equilibrium curve is presented. Some test functions 
are used in order to locate limit points and Hopf bifurcation points. 
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Two extended systems allow to calculate these points. The numerical 
results are presented in a bifurcation diagram. 

MSG: 37N25 37G10 37M20. 

keywords: limit point, Hopf bifurcation point, Luo-Rudy I model, arc- 
length-continuation method, Newton's method, computer program. 

1 Introduction 

The present paper is one of a series of research results, for the dynamical 
system associated to the Luo-Rudy I model, obtained under the coordination 
of Acad. Adelina Georgescu. 

Mathematical models of cardiomyocyte electrophysiology based on ex- 
perimentally determined kinetics of ion currents encompass over 40 years, 
since the early attempts of Denis Noble ([24], [25]) to accurate models of 
cell types in all regions of the heart that are now being incorporated into 
anatomically detailed models of the whole organ ([26j). These approaches 
introduced in successive steps several time-dependent ion current compo- 
nents, as well as a detailed dynamics of calcium in subcellular compart- 
ments, comprising release and reuptake from the sarcoplasmic reticulum 
([5j), specific calcium buffers ([16j), and the electrogenic Na/Ca exchanger. 
The Luo-Rudy I model of ventricular cardiomyocyte ( [22] ) , developed in the 
early 1990s starting from the Beeler-Reuter model ([T]), includes kinetics 
based on single-channel recordings. All current components are described 
by Hodgkin-Huxley type equations. This simplicity renders it adequate for 
mathematical analysis using methods of linear stability and bifurcation the- 
ory. 

Nowadays, there exist numerous software packages for the numerical study of 
finite - dimensional dynamical systems, for example MATCONT, CL_MATCONT, 
CL^MATCONTM ([3j, [l3]), AUTO [6j. In this paper, numerical results 
are obtained by using some new computer programs. 

2 Luo-Rudy I model 

In spite of its simplicity, which comes from the fact that it does not take 
into account earlier findings concerning Ca^^ dynamics, the Luo-Rudy I 
model [22] proved to be very realistic, incorporating data derived from 
single-channel recordings obtained during the 1980s with the advent of the 
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patch-clamp technique ( [15] , [23j ) . The model comprises only three time and 
voltage-dependent ion currents (fast sodium current, slow inward current, 
time-dependent potassium current) plus three background currents (time- 
independent and plateau potassium current, background current). The Luo- 
Rudy I model reproduces a wealth of experimental findings, like: the fast 
upstroke velocity of the action potential {Vmax = 400 V/sec), the behavior of 
the rising phase, late repolarization phase, and postrepolarization phase with 
changes in extracellular potassium concentration [K]o, monotonic Wencke- 
bach patterns and alternans at normal [K]o, nonmonotonic Wenckebach 
periodicities, aperiodic patterns, and enhanced supernormal excitability re- 
sulting in unstable responses and chaotic activity at low [K]o. 
Let us briefly describe how different experimental facts were taken into ac- 
count within the equations. The fast sodium current (Ing) incorporates both 
a slow process of recovery from inactivation and an adequate maximum con- 
ductance. The activation (m) and inactivation (h) rates are adapted from 
the Ebihara-Johnson Ijva model based on data from chicken embryo car- 
diac cells [7]. Two inactivation gates, fast and slow {h and j) were used to 
render it compatible with single-channel data proving that near threshold 
potentials sodium channels tend to open several times during a depolariza- 
tion (reopening phenomenon), and a significant fraction of channels do not 
open by the time of peak inward current. The start values of the slow in- 
activation gate j are obtained by setting joo = hoo, as suggested by Haas 
et al. [13] ■ The slow inward current (I si) is represented exactly as in the 
Beeler-Reuter model. The time-dependent potassium current (Ik) is con- 
trolled by a time-dependent activation gate {X) and a time-independent 
inactivation gate (Xj) with inward rectification properties, neither of which 
depends on [K]o, while the single-channel conductance is proportional to 
the square root of [K]o, as found in patch-clamp recordings on rabbit nodal 
cells |34j . The time- independent potassium current (Iki) is different from 
Iki of the Beeler-Reuter model, featuring two important properties dis- 
covered by Sakmann and Trube using patch-clamp methods ([2H], [22]): a 
square-root dependence of single-channel conductance on [K]o, and a high 
selectivity for potassium, as well as the inactivation gate Kl identified by 
Kurachi in single-channel experiments [19]. Since this current inactivates 
completely during depolarization, the model was supplemented with two 
other time-independent potassium current components: a [K]o -insensitive 
plateau current {Ixp), simulating the single-channel properties of the plateau 
current measured by Yue and Marban [36], and a background current (/(,) 
with a reversal potential Ei, = -59.87 mV. We should remark that, although 
derived from single-channel experiments, the conductances, gating and re- 
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versal potentials of these three current components were adjusted using a 
parameter estimation technique to fit the whole-ceU time-independent potas- 
sium current measured by Sakmann and Trube for different values of [K]o- 
The mathematical problem governing the membrane excitability of a ventric- 
ular cardiomyocyte, according to the Luo-Rudy I model ([22]), is a Cauchy 
problem 

u{0) = uo , (1) 
for the system of first order ordinary differential equations 

du , . 

-^ = Hv,u), (2) 
where u = (ui, . . . , ug) = {V, [Ca]i, h, j, m, d, f, X), r] = (t?i, . . . , 7713) 

= {1st, Cm, gNa, Qsi, QKp, Qb, [^(Ao, [Na\i, [K]q, [K]i, PRNaK, Eb, T), 

1 „ 

$1(7/,^) = [Ist + mU'iUAU^iui - ENa{'n7,'n8,m3)) 

+riiUQUj{ui - ci + C2 lnn2) 

+gK{ilw)Xi{ui){ui - -Ei<:(7?7,r/8,r?9, 7/10,7/11,7/13)^8 

+9X1 (??io) loo {m,mo,m'i,ui){ui - Eki (% , mo , ??i3 ) ) 

+r]^Kp{ui){ui - -Eii-p (779, 7/10, 7/13)) 7/6 (ui -7/12)], 

$2(??,'U) = -C2.rHUQUj{ui - Ci + C2lnti2) + C4(C5 - U2) , 

^i{r], u) = ai{ui) - {ai{ui) + Pi{ui))ui , £ = 3, . . . , 8 . 



m , 



The definitions of variables V, [Ca]i, h, j, m, d, /, X, parameters 1st, C, 
gNa, Qsi, QKi, QKp, Qb, [Na]Q, [iVa]i, [K]q, [K]i, PRNaK, Eh, T, constants 
ci, . . . ,C5, functions gx, E^a, Ek, Eki, Exp, Kl^o, Xi, Kp, ae, /3g, default 
values of parameters and initial values of variables in the Luo-Rudy I model 
are: V - transmembrane potential, [Ca]i - intracellular calcium concentra- 
tion, h and j - fast and slow inactivation variable of I^a (probability of gate 
h or j to be open), m - activation variable of Inu, d and / - activation and 
inactivation variable of Igi, X - activation variable of Ik, Xi - steady-state 
inactivation of Ik, K^oo - steady-state gating variable of Iki, Kp - steady- 
state gating variable of iKp, cti and - voltage dependence of opening and 
closing rates expressed as Boltzmann distribution functions for two distinct 
energy levels, 1st - steady depolarizing/hyperpolarizing applied current. Cm 
- membrane capacitance per unit area, g^a - maximal conductance of fast 
voltage-gated sodium current (per unit area), gsi - maximal conductance of 
slow inward (calcium) current, qk - maximal conductance of time-dependent 
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potassium current , gxi - maximal conductance of inward rectifier potas- 
sium current, gxp - maximal conductance of plateau potassium current, gi, 

- maximal conductance of background current, [A^a]o, [iVoJj, [K]q, [K]i - 
extra- and intracellular concentrations of sodium and potassium, PRNaK 

- sodium/potassium permeability ratio for Ik, Ej^a, Ek, Eki, E^p, Eb - 
reversal potentials of Ijya, Ik, Iki, Ixp, lb, T - absolute temperature. 

For the continuity of the model, the reader is referred to [21j . and for the 
treatment of the vector field $ singularities to [3j. $ is of class on the 
domain of interest. 

3 The one - parameter dynamical system associ- 
ated to the Luo - Rudy I model 

We performed the study of the dynamical system associated with the Cauchy 
problem ([I]), ([2]) by considering only the parameter r/i = 1st and fixing the 
rest of parameters. Denote \ = rji = 1st and 77* the vector of the fixed values 
of?72,...,r/i3. LetF:i2xM^M, F(A, u) = $(A, r?,, u), F = (Fi, . . . , Fg). 
Consider the dynamical system associated with the Cauchy problem ([T]), 
([3]), where 

'i=FiKu). (3) 

The equilibrium points of this problem are solutions of the equation 

F(A,u) = 0. (4) 

The existence of the solutions and the number were established by graphical 
representation in [3J, for the domain of interest. The equilibrium curve 
(the bifurcation diagram) was obtained in [3] , via an arc-length-continuation 
method ([llj) and Newton's method ([IQ]), starting from a solution obtained 
by solving a nonlinear least-squares problem ( [llJ ) for a value of A for which 
the system has one solution. In [3j, the results are obtained by reducing (j4]) 
to a system of two equations in {ui,U2) = iV, [Ca]i). Here, we used directly 

(SD- 

4 Arc-length-continuation method and Newton's 
method for (jlj) 

Let us write the arc-length-continuation method and Newton's method used 
to construct the equilibrium curve of (jl]). 
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Glowinski ([H], following H.B.Keller [17j . |18j ) chose a continuation equation 
written in our case as 

te + (^f = l. (5) 

1=1 

where s is the curvilinear abscissa. 

Let (A^=, It*) be a solution of dH obtained by solving a nonlinear least-squares 
problem ([H]) for a fixed value of A (A = A*) for which the system has 
one solution. To solve ([4]), let us consider the extended system formed 
by dH and ([5]), parameterized by s. Let As be an arc-length step and 
= u{nAs). We have the algorithm (following the case formulated in 
[TT]): take A(0) = A° = A*, n(0) = = and suppose that are 
given; for n > 0, assuming that A""""*^, u"^'^, A", are known, (X^^^ , u"^^^) 
is obtained by 

F(A'^+\n"+i) = (6) 

and 

huj -u^)^ + (A^ - A°)^^ = As if n = , 

E«^' - <) ^ + (A"+^ - A") ^ (7) 

= As if n > 1 . 



In order to calculate '^^j^^ , we obtain the following relations. FromdH), 

we have 

^ 9F,(A°,^°) dn,(0) , dFj{XO,n^) dX{0) 

+ — dx dT-"'^-^'---'^- 



Let 



a!ni(0) „ (iA(O) . ^ . ... 
-Ui— — ,z = l,...,8. (9) 



ds ds 
u = (ni, . . . , -Us) is the solution of 

.8. ,10) 
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Prom ([5]), we have 

ds 



E^^ + i)(^)' = i- (11) 



i=l 

In dZl), let us denote u* = u^, X* = A°, u** = X** = ^ if n = and 

u* = u", A* = A'^, u** = ""^"J"' , A** = if n > 1. 

The algorithm which we use to construct the branch of solutions for is 
the following: 

1. given A*^, u^, solve ([TO]) to obtain u; 

2. obtain ^ from (HH) and ^ from ([9]); 

3. A(0) = A*^ = A* and u{0) = = are taken as discussed above; for 
n > 0, taking A", u" as initial iteration, the following algorithm based on 
Newton's method calculates A""*""*^, n""^^ : obtain (A™'"'""^, ■u'""'"^) using 

^ aF,(A",n"-) ^^+1 , aF,(A"^,0 
^ 9^1^ + aA ^ 

2=1 * 



_ ^dFj{X\vI^ „ 9Fj(A'",u^ 

~ hi ^ " 

-F,(A-,0,J = l,---,8 (12) 

8 8 

^ -urnf +1 + A**A™+i = u*ur + A* A** + As ; 

1=1 i=l 

calculate the eigenvalues of the Jacobian matrix DuF{X'^~^^ , u"^~^^) by the QR 
algorithm, calculate ■0lp(A"+-^, ti"+-'^) by (fT3|) . and calculate ■ipni^"''^^ ,u"'~^^) 
by (HID; 

4. the algorithm is stopped after an imposed number of iterations for n. 



5 Limit points and Hopf bifurcation points 

In order to locate the limit points and the Hopf bifurcation points on the 
equilibrium curve of ([1]), two test functions ([20], [E], [32], [33]), tpip and 
tpH, are evaluated at each iteration (A""*"^, ■u""'"^) of the algorithm from the 
end of section IU where 

2PLp{X,u) = det{D^F{X,u)), (13) 



iPh{X,u) = det{2DuF{X,u) Ig) 



(14) 



8 



Bichir, Georgescu, Amuzescu, Nistor et al. 



ipH is evaluated using formula (jlSp . For a, n x n matrix A with elements 
{ttij}, the following m x m matrix, m = ^n{n — 1), is obtained ([2Q], [12]) 
based on the definition of the bialternate product of two matrices, 



(2^0 In){p,q),{r,s) 



^ps 1 



^pp 



if r = ^ , 
if r 7^ p and s = q , 
+ aqq , if r = p and s 



'•qr 1 



^qq , 

if r = p and s ^ q . 
if s=p, 



(15) 



, otherwise . 



The rows are labeled by the multi-index {p,q) (p = 2,3, ... ,n, q = 1,2, .. . ,p— 
1), and the columns are labeled by the multi-index (r, s) (r = 2,3, ... ,n, 
s= l,2,...,r-l). 

If V'LP has opposite signs at two points (A",ti") and (A"+^,n"+^), then a 
limit point exists between (A",n") and {X^'^^ ,u"'^^). If TpH has opposite 
signs at these two points, then it is possible that a Hopf bifurcation points 
exists between them. The existence of a Hopf bifurcation point is decided by 
studying the form of eigenvalues of DuF{X^,u'^) and of DuF (X^^^ , u^^^) , 
since ipn can be zero if there is a pair of real eigenvalues of opposite sign 
and with equal modulus. (We had in view the existence of a pair of complex 
conjugate eigenvalues, for each matrix, and a change of the sign of their real 
part). 

For the cases where the test functions detect a limit point or a Hopf bi- 
furcation point, we retain the results of one of the iterations (A",n") or 
(A""^^, li""''^) of the algorithm presented at the end of section [H namely the 
iteration where the modulus of the test function is smaller. These results 
are an initial iteration for Newton's method applied to one of the equations 

[32], m) 

G{X,u,h)=0, (16) 
H{\,f3,u,h,g)=0. (17) 



The components {X,u) of the solution of (I16p represent a limit point of (|4]). 
The components (A, u) of the solution of (fT7|) represent a Hopf bifurcation 
point of H and G are defined by 



F{X,u) 

DuF{X,u)h 

hk-l 



(18) 
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H : R^-^^"^ H{\, /3, u, h, g) 



F{X,u) 

DuF{\u)h + pg 
DuF{X,u)g - l3h 
hk-l 
9k 



(19) 



where /c is a fixed index, 1 < A; < 8. The extended system (|17|) determines 
a Hopf bifurcation point {X,u), a pair of purely imaginary eigenvalues ±/?i 
of the Jacobian matrix in (A, u), and a nonzero complex vector h + ig. 
Newton's method applied to the equation (fT6|) is: let = (A*^, u^, h^) be an 
initial iteration, where A", are as specified above and h? = (1,0, 0, 0, 0, 0, 0, 0); 



= 1; for m > 0, calculate v'^+^ = (A'^+\ 



u 



m+1 ^m+1 



usmg 



DG{v"'){v 



m+1 



(20) 



Newton's method applied to the equation (fTTj) leads to: let = (A*^, , 

g^) be an initial iteration, where A*^, are as specified above, 
is the positive complex part of one of the two complex conjugate eigen- 
values of the Jacobian matrix in (A'^, n"), = (1,0,0,0,0,0,0,0), g^ = 
(0,0,0,0,0,0,0,0); = 1; for m > 0, calculate w^+i = (A'"+\ /3'"+\ u^+S 
/jm+i^ 5™+i) using 

DH{w"'){w"'+^ - y;"^) = - Hiw"^) . (21) 



6 Numerical results 

Based on the computer program for [2], the algorithm at the end of sec- 
tion H] was transformed into a new computer program. Two new computer 
programs were written in order to solve (jl6p and (|17p by (j20p and ()2ip re- 
spectively. 

We took Cm = 1, 9Na = 23, gsi = 0.09, gx = 0.282, gKi = 0.6047, 
gxp = 0.0183, Gb = 0.03921, [Na]o = 140, [Na]i = 18, [i^]o = 5.4, 
[K]i = 145, PiiATaX = 0.01833, Eh = -59.87, T = 310. 

The equilibrium curves (the bifurcation diagram) are presented in figure [T] 
for the domain of interest. The variation of the functions tpLP and ipn Sire 
also represented. "LPl" and "LP2" indicate two limit points. "H" indicates 
a Hopf bifurcation point. 




Figure 1: The bifurcation diagram and the variation oi ipLP and tpn- There 
exist two Umit points (LPl, LP2) and a Hopf bifurcation point (H). 
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The behavior of the dynamical system changes with the value of the variable 
parameter Igf The two limit points separate three branches of stationary 
solutions (first graph in Fig. 1). While solutions on the middle branch (ac- 
cording to values of V) are always unstable and those on the lower branch 
are always stable, on the upper branch there is a region where the system 
features oscillatory behavior. Oscillations are either damped, at the left of 
the Hopf bifurcation point, or amplified until the system falls on the lower 
branch of solutions to the right of the Hopf bifurcation. Amplified oscil- 
lations represent early afterdepolarizations (EADs), a condition prone to 
result in life-threatening arrhythmias. A recent study of the Luo-Rudy I dy- 
namical system for variable relaxation time constants of the gating variables 
d, /, and X ([35J), has proved that oscillations resulting in EADs appear 
above a Hopf bifurcation point for a fast subsystem, comprising the variables 
V, d, and /. Moreover, these EADs can result in chaotic behavior when the 
system is paced at a constant cycle length. The same group has shown on 
detailed three-dimensional ventricular electrophysiology models that EADs 
occurring in certain regions can synchronize, resulting in polymorphic ven- 
tricular tachycardia or torsades-de-pointes ([SOj). 

In conclusion, our study, focused on analysis of the Luo-Rudy I system as a 
whole in conditions of variable parameter 1st, has identified two limit points 
and a Hopf bifurcation point, separating different regions of stability, some of 
them featuring amplified self-sustained oscillations defined as EADs on the 
time trajectories, and which may result in dangerous ventricular arrhythmias 
by synchronization. In contrast to the majority of previous arrhythmoge- 
nesis studies, which attributed the generation of this phenomenon to an 
individual condition, such as altered gating of an ion channel type due to 
gene mutations or modulation by physiological or pharmacological mecha- 
nisms, our results prove that arrhythmogenesis may result as an emergent 
feature of the system as a whole, and not of its individual components. 
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